Chiral Structure of F-actin Bundle Formed by Multivalent Counterions 
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The mechanism of multivalent counterion-induced bundle formation by filamentous actin (F-actin) 
is studied using a coarse-grained model and molecular dynamics simulation. Real diameter size, 
helically ordered charge distribution and twist rigidity of F-actin are taken into account in our model. 
The attraction between parallel F-actins induced by multivalent counterions is studied in detail and 
it is found that the maximum attraction occurs between their closest charged domains. The model 
F-actins aggregate due to the like-charge attraction and form closely packed bundles. Counterions 
are mostly distributed in the narrowest gaps between neighboring F-actins inside the bundles and the 
channels between three adjacent F-actins correspond to low density of the counterions. Density of the 
counterions varies periodically with a wave length comparable to the separation between consecutive 
G-actin monomers along the actin polymers. Long-lived defects in the hexagonal order of F-actins in 
the bundles are observed that their number increases with increasing the bundles size. Combination 
of electrostatic interactions and twist rigidity has been found not to change the symmetry of F- 
actin helical conformation from the native ^ symmetry. Calculation of zero-temperature energy 
of hexagonally ordered model F-actins with the charge of the counterions distributed as columns 
of charge domains representing counterion charge density waves has shown that helical symmetries 
commensurate with the hexagonal lattice correspond to local minima of the energy of the system. 
The global minimum of energy corresponds to symmetry with the columns of charge domains 
arranged in the narrowest gaps between the neighboring F-actins. 
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I. INTRODUCTION 

F-actin is a highly charged polyelectrolyte that plays a 
key role in the cytoskeleton of eukaryotic cells, where it 
forms a network of bundles made by crosslinking proteins 
It has also been observed that multivalent counteri- 
ons can cause these macromolecules to attract each other 
(the phenomenon known as like-charge attraction) and 
form aggregates 0-3 . 

When the separation between two polyelectrolytes like 
F-actin is considerably larger than their diameter, they 
can be modeled as uniformly charged lines, cylinders, or 
bead-spring chains, for the calculation of their electro- 
static interaction potential [5|-|25|. Despite their sim- 
plicity, such coarse-grained minimal models have been 
shown to suffice for obtaining the essence of the overall 
behavior of these polyelectrolytes in solution. Analyti- 
cal and computational studies using these models have 
shown that the interaction potential between two simi- 
larly charged stiff polyelectrolytes in multivalent salt so- 
lution crucially depends on their separation, mutual ori- 
entation and the concentration and valence of the salt 
[H, HO, Il6'-'28'l. Many aspects of the collective behavior 
of stiff polyelectrolytes in multivalent salt solution, such 
as the formation of aggregates of various structures and 
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high sensitivity of these structures to the salt concen- 
tration and valence, are believed to originate from the 
complexity and high anisotropy of their interaction po- 
tential 18, 20, 27, 28]. The complexity is the result of the 
combination of the highly anisotropic shape of stiff poly- 
electrolytes and the simultaneous presence of short- and 
long-range interactions in the system, and is also known 
to lead to frustration and slow dynamics [H, Is^l • 

The interaction between these macromolecules at close 
separations, however, cannot be studied using such sim- 
plistic models. In this case, the detailed structure of the 
macromolecules such as their helical symmetry, charge 
distribution, thickness, dielectric constant, twist rigidity 
and so on, will need to be taken into consideration. It 
has been observed that like-charge attraction between he- 
lical polyelectrolytes such as F-actin and DNA induced 
by multivalent salt could organize them in closely packed, 
hexagonally ordered bundles of parallel polyelectrolytes. 
Native helical symmetry of these polyelectrolytes ma y be 
incommensurate with the mentioned hexagonal order |3l|- 
[S^l- In this case, the interplay between twist rigidity 
of the polyelectrolytes and electrostatic correlations de- 
termines their equilibrium conformation in the bundle. 
There have been a number of theoretical studies on in- 
teraction between helical macromolecules in the litera- 
ture [lli [35|. Experimental study of bundle formation 
by F-actin in the presence of multivalent salt has shown 
that there are major differences between F-actin bundles 
formed in this way and those formed non-electrostatically 
[1]. To study like-charge attraction between such poly- 
electrolytes in close separations the model should contain 



2 



details of the charge distribution on the polyelectrolytes 
and their orientational and torsional degrees of freedom. 
Coarse-grained models (such as the bead-spring model) 
are not suitable for such studies. On the other hand, 
computer simulation of these macromolecules in atom- 
istic scale is too time consuming, while a large number 
of the internal degrees of freedom will also be redun- 
dant and have no role in the interaction between the 
macromolecules. Therefore, a sufficiently detailed coarse- 
grained model that incorporates the helical structure of 
the macromolecules seems to be the most appropriate 
starting point for such a study. 

In this paper, we study the mechanism of like-charge 
attraction between coarse-grained model F-actins in the 
presence of multivalent counterions using molecular dy- 
namics (MD) simulations. We find for two parallel 
F-actins that the maximum attraction occurs between 
monomers whose charged domains are in the closest sep- 
aration. The attractive potential of mean force between 
F-actins reveals no detectable dependence on their helical 
symmetry. The like-charge attraction between F-actins 
organizes them in bundles of hexagonal lattice structure 
containing long-lived defects. The number of the defects 
in a bundle is found to be an increasing function of its 
size, namely the number of F-actins forming the bundle. 
In the cross section of each bundle, the low- and the high- 
counterion-density regions are observed that correspond 
to the narrowest gaps between neighboring F-actins and 
the centers of triangles formed by centers of three adja- 
cent F-actins, respectively. Competition between electro- 
static interactions and twist rigidity of F-actins appears 
to be unable to change their native helical conformation 
of ^ symmetry into structures with other symmetries 
(the definition of F-actin helical symmetry is given be- 
low, see Sec. |TT|. The helices of the F-actins are observed 
to form domains of uniform phase that are connected via 
defects analogous to those in 2D XY model. Consider- 
ing the charge of counterions distributed as columns of 
frozen charge domains similar to the counterion charge 
density waves introduced in Ref. Q, we calculate zero- 
temperature energy of the system as a function of the 
helical symmetry of the F-actins. We find that there 
are numerous local minima in the energy of the system 
corresponding to helical symmetries commensurate with 
the hexagonal order of the F-actins in the bundle. The 
ground state of the system corresponds to the |y heli- 
cal symmetry of F-actins when the columns of charged 
regions are at the midpoints between the axes of neigh- 
boring F-actins. 

The rest of the paper is organized as follows. The 
model and the simulation method are introduced in Sec. 
im The study of like-charge attraction between a pair of 
F-actins is presented in Sec. IIIII In Sec. IIVI the counte- 
rion density profile, as well as the positional and helical 
structures of the filaments in the bundle are studied. The 
helical symmetry of F-actins in the bundle is studied in 
Sec. |V] using the simplified smeared counterion distribu- 
tion model. Section IVll concludes the paper. 



II. THE MODEL AND THE SIMULATION 
METHOD 

We use a coarse-grained model for F-actin that takes 
into consideration a number of key macromolecular fea- 
tures such as its diameter (that is relatively large, say as 
compared to the mobile ions) , helical order of its charge 
distribution, persistent length, rotational and twist de- 
grees of freedom, and its twist rigidity. We construct the 
model F-actin by assembling two kinds of spheres (see 
Fig. [T]): a large neutral sphere that makes the back- 
bone of F-actin (that represents the volume of actin) 
and a smaller negatively charged sphere that carry the 
net charge of F-actin (which has a linear charge den- 
sity -actin — — ^). Each G-actiu monomer is thus 
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modeled as a set of a large sphere and an attached small 
charged sphere. The adjacent large spheres are bonded to 
each other by a harmonic potential, = ^kh{r — vq)^ , 
of spring constant kh — 500cr~^, where tq is the aver- 
age separation between them and a is the MD length 
scale (see below) . Three-particle angle potential, = 
kbendi^ — COS0), with khend = 1000 is uscd to mimic F- 
actin bending rigidity and keep our model F-actins as 
rod-like polymers (since their lengths are considerably 
smaller than the persistence length of F-actin). The po- 
sitions of the small spheres are determined from X-ray 
observation-based four-sphere model of G-actin |36l - [38| . 
These spheres play the role of subdomain-1 (sdl), which 
is known to carry most of the G-actin charge. According 
to the four-sphere model, the position of the sdl domain 
of the nth monomer in the regular conformation of F- 
actin, can be given by 

p = 28.331 

0n = (f>l + {n~ 1)00 (1) 
Zn = zi + {n- l)zo, 

in cylindrical coordinate system. Here, zi and 4>i are the 
height and the azimuthal coordinates of sdl of the first 
monomer and 0o a-nd zq are the increments of 4> and z 
coordinates from each monomer to the next one. We set 
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FIG. 1: The model F-actin constructed by two kinds of 
spheres. The large overlapping blue spheres model F-actin 
backbone and carry no charge. The yellow spheres carry all 
the net negative charge of G-actin monomers. The small red 
sphere (whose diameters are scaled by a factor of three to 
make them visible) represent the counterions. The spatial or- 
der of the yellow spheres around the backbone is determined 
from the real helical order of F-actin and the structure of 
G-actin monomer obtained from X-ray experiments pj6-,38.] . 
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zq = 28.7 A and determine the value of (po from the heh- 
cal symmetry of F-actin, such that for heUcal symmetry 
of — , 00 — — "^^^ ■ For example, for helical symmetry of 
fUo = -166.2°. Small spheres in the model 

F-actin are bonded to the large ones by a harmonic po- 
tential of spring constant fc/j. To model twist rigidity 
of F-actin ^ (ntwist = 8 x lO'^^TVm^), the dihedral 
potential = ktwist [1 — cos(0 — (j)o)] is applied be- 

tween two neighboring monomers with kt^ist — 6900, 
where (p — (f)^ is deviation of dihedral angle between two 
consecutive monomers from its native value, (pQ. From 
the known charge density of F-actin, it can be calculated 
that each small sphere carries the charge ofq = — lie (eis 
the elementary charge) . The charge of each small sphere 
is considered as a point charge in its center. We explicitly 
add 3-valent counterions of charge 3e as Lennard-Jones 
(LJ) particles of diameter a = 2.2A (MD length scale 
in our simulations) to neutralize the charge of F-actins. 
Lennard-Joncs potential also introduces an energy scale, 
e, to the system. We use 3-valent (not 2-valent) counteri- 
ons to strengthen the electrostatic correlations in the sys- 
tem and probe the role of these correlations easily. One 
should note that because of the coarse-graining method 
in our model there is no essential difference between 2- 
and 3-valent counterions. 



Considering that in typical actin bundling experiments 
3 the salt concentration is 36-108 niM, electrostatic in- 
teractions are strongly screened and act only at short 
ranges (the screening length is Aocbyc — lOA). Like- 
charge attraction between polyelectrolytes is known to 
originate from correlations and fluctuations of multiva- 
lent counterions in their vicinity. To account for the ef- 
fect of salt and avoid the time consuming calculations of 
the long-ranged Coulomb interactions, we use a screened 
Debye-Hiickel potential and explicitly add counterions to 
the system. The Debye-Hiickel potential for the electro- 
static interaction between the charged particles i and j is 
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where I b — 7 A is the Bjer- 



rum length (in water at room temperature), Zi(Zj) is the 
valence of charged particle and rij is the separation 
between charges i and j. We use = lOcr. One should 
note that effects such as electrostriction which may have 
some contribution [40l - l42j are not considered here. We 
denote the numbers of F-actins, 3-valent counterions and 
monomers of each F-actin by Np, Nc, and Nm, respec- 
tively. Our MD simulations are performed with the MD 
simulation package ESPResSo ^43]. Periodic boundary 
condition (PBC) is applied and a Langevin thermostat 
with friction coefficient F = 1.0-^/^ is used to keep the 
temperature fixed at fc^T — e (m is the mass of the 
mobile ions). The simulation box is rectangular and its 
length in z direction is equal to the length of F-actins 
to minimize the end effects. The lengths of the box in 
the other two directions are the same and determined 
by considering the cut-off length of the screened elec- 
trostatic interaction and the number of F-actins in the 
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FIG. 2: The potential of mean force per monomer between a 
pair of parallel F-actins of the native helical symmetry and 
two other symmetries, || and in the presence of 3-valent 
counterions as a function of their separation, R. Each data 
point is obtained from averaging over 10^ integration time 
steps after equilibration of the system. Error bars which are 
shown only on data points of || symmetry are obtained from 
averaging over 6 pairs (central F-actin with 6 others surround- 
ing it). As it can be seen, the values of the potential of mean 
force for the three helical symmetries are the same within the 
error bars. Inset: Top view of the configuration of 7 F-actins 
and average forces acting on their monomers for a given value 
of R. 



system. Truncated and shifted LJ potential. 
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if r < To + Tc 

if r > To -|- 
(2) 

is used to model excluded volume interactions. Tj. = 
and To are the cut-off and offset of this potential |43{, 
respectively. In the cases of backbone and sdl spheres 
of model F-actins, the corresponding diameters of the 
spheres, namely df, = 28.33A x 2 = 56.66A and asdi = 
I6.8IA X 2 = 33.62A, are used to set the value of the 
variable rg in the LJ potential. In the case of counterion- 
counterion excluded volume interactions, ro = 0. LJ time 

scale in our simulations is tq — and the integration 

time step is r = 0.02to. 



III. LIKE-CHARGE ATTRACTION BETWEEN 
TWO PARALLEL F-ACTINS 

To calculate the potential of mean force between a 
pair of parallel model F-actins in the presence of 3-valent 
counterions, we first fix Np — 7 of them parallel to each 
other on a hexagonal lattice of separation R as shown 
in the inset of Fig. [51 Keeping all of their degrees of 
freedom blocked, we leave the counterions to equilibrate. 
We then release the rotation and twist degrees of freedom 
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{(j>i and (j>Q in Eq. [ij of all F-actins, such that they can 
explore their equilibrium conformations. After equilibra- 
tion of the system, we calculate the time average of the 
interaction forces between the central F-actin and 6 oth- 
ers surrounding it. From integration of the mean force 
(averaged over time and over 6 pairs of F-actins) with 
respect to the separation R, we calculate the potential of 
mean force for a pair of F-actins. The potential of mean 
force per monomer versus the separation is shown in Fig. 
m We calculate the potential of mean force for three dif- 
ferent values of ipQ corresponding to helical symmetries 
and 11 . Our results show that the interaction 
potential as a function of the separation does not depend 
on F-actins helical symmetry and no detectable differ- 
ence for above-mentioned symmetries is observed (see 
Fig. [5]). We also observe that the difference between 
initial phases, 0i, of the interacting F-actins tends to 
vanish in equilibrium. 

To understand the mechanism of interaction between 
a pair of neighboring model F-actins it is instructive to 
look at the forces acting on each of their monomers. The 
average forces acting on monomers of two interacting F- 
actins are shown in Fig. |31 which shows that the mu- 
tual attraction is strongest between the few monomers in 
each helical period whose charged domains are closest to 
each other. Interactions between the other monomers are 
found to be considerably smaller than that of the over- 
lapping charged domains. The few overlapping domains 
produce a strong attractive potential for the counterions 
and gather them in the narrow space between these do- 
mains. 
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FIG. 3: The time average of the forces acting on monomers of 
two parallel model F-actins in the presence of 3-valent coun- 
terions. Equilibrium configuration of the F-actins is shown in 
the upper panel. The strongest attractions are between the 
monomers whose charged domains are closest to each other. 
The values of the variable for the two F-actins are the 



IV. BUNDLE FORMATION BY F-ACTINS AND 
ITS MECHANISM 

To study the bundle formation of the filaments, we 
first fix Np parallel F-actins of length Nm = 26 on a 
hexagonal lattice of spacing R = 80 A°, larger than the 
value corresponding to the minimum of the potential of 
mean force (see Fig. [2]). We then leave the counterions 
to equilibrate, and subsequently release all the degrees 
of freedom of the F-actins, except that we constrain the 
two ends of the filament to move on the two boundary 
surfaces of the simulation box. We have also repeated 
the same procedure with non-hexagonal initial configu- 
rations of F-actins. We observe that regardless of the ini- 
tial configuration of the F-actins they attract each other 
and form a bundle in which their axes form a hexagonal 
lattice. Moreover, the initial phases of the helices (vari- 
able (f>i in Eq. [T]) are observed to approach each other in 
the course of the aggregation dynamics. 

To visualize this behavior we represent the initial 
phase, (j>i, of each F-actin by a vector directed from F- 
actin axis to sdl of its first monomer and probe the direc- 
tional ordering of these vectors in the x — y plane in our 
simulations. We perform these simulations for different 
system sizes, Np — 7 — 52. In the case of Np = 7, the 
system easily finds its low-energy state in which parallel 
F-actins are ordered hexagonally and their phase vectors 
are mostly aligned in the same direction. For large Np, 
however, usually there are some defects in the hexago- 
nal order of F-actins in the bundle, with their number 
increasing when the number of F-actins in the system 
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FIG. 4: a) Initial configuration of A'p = 30 parallel F-actins 
(top view) with non-ordered phase vectors, b-f) Five typical 
configurations of F-actins after running the system for a long 
time (~ 5x 10^ integration time steps), b) Both positional and 
directional orders are approximately complete, c-f) Different 
kinds of defects in F-actins positional hexagonal order and 
directional order of their phase vectors. 

increases. These defects are very long-lived and for the 
systems larger than Np ~ 10 it is not possible to observe 
their disappearance and to find defect-less hexagonal or- 
der in the system in our simulations. Moreover, defects 
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FIG. 5: The order parameter, X, defined in Eq. [S] as a 
function of time obtained from simulation of A'p = 30 model 
F-actins. Through directional ordering of the phase vectors 
the order parameter changes from zero to a saturation near 
unity. Typical snapshots of the phase vectors are shown for 
randomly oriented and orientationally ordered configurations. 

in the long-range directional order of the phase vectors 
(similar to those found in 2D XY model) are also ob- 
served in the system (see Fig. S]). The positional and 
the orientational defects are strongly correlated. Namely, 
when there are positional defects in the system, certainly 
orientational defects are also present. Orientational de- 
fects are also observed in systems without any positional 
defects. Low-energy state of the system corresponds to 
a configuration in which F-actins form a hexagonally or- 
dered bundle with no defect and with their phase vectors 
all aligned in the same direction. We define an order 
parameter, X, as 

X = S<,,>cos(01-0l), (3) 

in which the summation runs over nearest neighbor pairs 
to probe directional ordering of the phase vectors. A 
typical behavior of the order parameter, X, as a function 
of time is shown in Fig. [S] The figure shows that the 
order parameter starts from the initial value of zero cor- 
responding to complete lack of order in the phase vectors 
orientation, and saturates to a value that represents the 
degree of alignment between the phase vectors. 

Simulations of F-actins with the value of 0o set to the 
native ^ symmetry show no appreciable change in the 
F-actins helical symmetry upon bundle formation. Al- 
though local deviations of dihedral angles between adja- 
cent monomers along F-actins are observed, the average 
value of the dihedral angles agrees well with the native 
value. We performed these simulations with two other he- 
lical symmetries of the model F-actins, namely H and , 
and observed no noticeable changes in the results. Statis- 
tics and averages of deviations of dihedral angles between 
adjacent monomers from their native values for three dif- 
ferent helical symmetries are the same within the error 
bars. We also performed simulations with model F-actins 
with native helical symmetry of ^ but using initial con- 
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FIG. 6: Top view of counterion density profile inside a hexag- 
onally ordered bundle of A'p — 30 model F-actins (a part of 
the system is shown here). The counterion density has its 
maximum value in the middle of the separation between the 
axes of neighboring F-actins. The centers of triangles formed 
by three neighboring F-actins have the lowest value of the 
counterion density. 



formations (after equilibration of the counterions) that 
had symmetries of ^ and (two separate simulations). 
We found that in both cases, after releasing the F-actins 
degrees of freedom their symmetry returned back (in a 
short time) to the native These observations suggest 
that electrostatic correlations are not capable of keeping 
changing the helical symmetry of F-actins in bundle. 

To study the counterion distribution and the force field 
acting on G-actin monomers in a bundle of hexagonally 
ordered parallel F-actins, we construct a bundle of Np 
fixed F-actins with their axes on a hexagonal lattice and 
leave the counterions to equilibrate. We choose the lat- 
tice constant as the value that corresponds to the mini- 
mum of the potential of mean force shown in Fig. [2] We 
then release the twist and rotational degrees of freedom 
of F-actins and allow them to equilibrate. After equili- 
bration of the system we probe the counterion distribu- 
tion and the forces acting on the monomers of F-actins 
and any deviation of the twist order of F-actins from its 
native symmetry. A typical counterion density profile 
chosen from a part of the system containing Np = 30 
model F-actins is shown in Fig. [Bl in which the distri- 
bution profile is obtained using the superposition of 10'^ 
snapshots of the system chosen from 10^ integration time 
steps of the simulation. 

We also calculate counterion-counterion pair correla- 
tion function in z direction. This correlation function 
shows that the most repeated separation between be- 
tween counterions in z direction is the separation between 
sdls of consecutive monomers along F-actins (see Fig. 
[7]). This result combined with the result shown in Fig. [S] 
shows that the counterions are mostly distributed in the 
narrowest gap between two neighboring F-actins inside 
the bundle as domains of separation equal to monomer- 
monomer separation along the F-actins. The centers of 
the channels formed by three adjacent F-actins corre- 
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FIG. 8: Forces acting on monomers of A'p — 7 model F-actins 
in the presence of 3-valent counterions (top view) averaged 
over 10'^ snapshots of the system chosen from 10^ integration 
time steps for hehcal symmetries ^ (a) and ^ (b). 



FIG. 7: Pair correlation function of the counterions along z 
axis inside a hexagonally ordered bundle of parallel F-actins 
which its cross-section in x ~ y plane is shown in Fig. (6] 
As it can be seen the most repeated separation between the 
counterions along the F-actins is approximately equal to the 
separation between consecutive G-actin monomers. 

spond to a lower density of the counterions despite the 
picture of charge density wave presented in Ref. Q. 

Time-average forces acting on the actin monomers in 
a bundle of A'p = 7 model F-actins for two values of 
00 corresponding to and symmetries are shown in 
Fig. [5] The main difference between helical symmetries 
^ and YY is related to their commensurability with the 
hexagonal lattice. Despite the native ^ symmetry of F- 
actin, symmetries like H ff commensurate with 
hexagonal lattice (considering its definition, symmetry ^ 
is commensurate with hexagonal lattice if m is divisible 
by 6). As Fig. [S] shows, the strengths and directions of 
the forces that act on actin monomers are slightly more 
regular in the case of . Apart from that, we observed 
no detectable difference between the time-average forces 
acting on the F-actins for the three different symmetries. 
We note that we have found some evidence suggesting 
that the number (or density) of positional defects for the 
same value of A'p is higher for the native symmetry of 
as compared to those that are commensurate with the 
hexagonal order of F-actins. Although this observation 
seems plausible, it needs to be probed more systemati- 
cally using larger systems. 



V. HEXAGONALLY ORDERED F-ACTIN 
BUNDLE WITH COLUMNS OF CHARGE 
DOMAINS OF OPPOSITE SIGN 

Here we consider a bundle of hexagonally ordered par- 
allel F-actins with the neutralizing charge of the counte- 
rions distributed as equally charged domains inside the 
bundle. The charge domains form one dimensional lat- 
tices of spacing 2 x zq = 57 AA along F-actins and are 
regularly distributed inside the bundle. Considering sym- 



metries of the system, two different arrangements of the 
columns of charged regions are possible. Top view of 
these two arrangements are shown in Fig. [9] The amount 
of the charge of these domains, which are considered as 
point charges in centers of gray spheres (see Fig. |9|), are 
determined in a way that the total charge of the sys- 
tem is zero. For both of these arrangements we calculate 
total energy of the system composed of electrostatic in- 
teractions (Debye-Hiickel potential with Xoebye = 22A) 
and twist energy of F-actins, for different values of regu- 
lar deviation of F-actins helical order from the native ip- 
6 

symmetry. In Fig. I10| energy of the system versus the 
value of the deviation from the native symmetry for the 
separation between neighboring F-actins, R — 60A, is 
shown. As it can be seen, there are some local minima in 
the energy of the system corresponding to symmetries of 
F-actin which are commensurate with hexagonal lattice. 
The global minimum corresponds to symmetry of ar- 
rangement (b) in Fig. |9l Also, as it can be seen in this fig- 
ure, energy of the system is considerably lower and local 
minima are considerably deeper in the case of arrange- 
ment (b) relative to arrangement (a). A similar study for 
arrangement of multivalent counterions inside a bundle of 
polyelectrolytes modeled as uniformly charged lines has 
been performed [2^. Despite the model studied in Ref. 
[25} . helical order of F-actin charge distribution enters 
the subject of commensurability into the problem in our 
model. We also looked at the effect of the slide of the 
neighboring one-dimensional lattices of charge domains 
with respect to one another and found that such slides 
increase the system energy and the low energy state of 
the system corresponds to the configuration with no slide 
of these columns relative to one another. Considering 
these results in addition to the counterion distribution 
obtained from the bulk simulations at room temperature 
(Fig. [5]), it seems that the counterions have a high ten- 
dency to the narrowest gap between two F-actins in the 
bundle and the picture shown in Fig. |6]is the dominant 
scenario relative to the other picture in which counteri- 
ons are mostly distributed inside the channels formed by 
three adjacent F-actins. 
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(b) 
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FIG. 9: Top view of two possible arrangements of the columns 
of charge domains inside the F-actins bundle introduced in 
the text. Counterion charge domains and sdls of G-actin 
monomers along the F-actins are shown by gray and black 
circles, respectively. 
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VI. CONCLUSIONS AND DISCUSSION 

In conclusion, using a coarse-grained model and MD 
simulations, the mechanism of like-charge attraction be- 
tween F-actins and bundle formation in the presence of 
multivalent counterions have been studied. Important 
macromolecular features of F-actin such as its diameter 
size, persistence length, helically ordered charge distri- 
bution, and twist rigidity have been considered in our 
model. It has been shown that parallel F-actins attract 
each other in the presence of 3-valent counterions and 
the maximum attraction occurs between monomers that 
their charged domains are in the closest separation. Due 
to the like-charge attraction, F-actins aggregate and form 
closely packed bundles of parallel F-actins. It has been 
found that the phase difference between helical structure 
of neighboring F-actins in a bundle tends to vanish. De- 
fects in the long-range orientational order of the phase 
vectors with interesting similarities to those of 2D XY 
model have also been observed. In addition, long-lived 
positional defects in the hexagonal order of F-actins in 
the bundles have been observed that their number in- 
creases with increasing the bundle size. Existence of 
similar defects have also been reported in experimental 
studies of similar systems 0, |^ . 

Counterions inside the bundles have been observed 
that are mostly distributed in the narrowest gap between 
neighboring F-actins and the channels formed by three 
adjacent F-actins in the bundles contain lower density of 
the counterions. The density of counterions in the nar- 
rowest gap between neighboring F-actins in the bundles 
varies periodically along F-actins with a wave length ap- 
proximately equal to the separation between consecutive 
G-actin monomers along the actin polymer. 

Our results obtained from simulations presented in Sec. 
IIVI show that in the bundle formation process, deviation 
of dihedral angle between consecutive monomers along 
F-actins from its native value happens. It seems that by 
these local deviations the system decreases the incom- 
mensurability effects of F-actins native structure with 
the hexagonal order. Such local deviations from a regu- 
lar symmetry in F-actin helical structure has also been 



FIG. 10: Energy per monomer of the central F-actin in ar- 
rangements (a) and (b) of Fig. [9] versus the deviation of F- 
actins helical symmetry from the native ^ symmetry, (j) — 4>o- 
Minima in the energy curve correspond to symmetries of F- 
actin which are commensurate with hexagonal order. The 
global minimum corresponds to symmetry . F-actin in this 
symmetry is 1.2° undertwisted relative to its native symme- 
try. Tr refers to room temperature. 



reported in experimental studies [4J]. Changing of F- 
actin helical symmetry from ^ to other symmetries has 
not been observed in our simulations. Positional defects 
of life time longer than our simulation times have been 
observed in the hexagonal order of F-actins in the bun- 
dles. The number of these defects has been found that 
increases with increasing the bundle size. Also, it has 
been observed that no considerable coarsening of these 
defects happens in our simulations because of very slow 
dynamics of their evolution inside the bundle. Although 
additional studies of these defects in larger systems are 
needed, existence of such defects could be a reason for fi- 
nite sided bundles of F-actin in multivalent salt solution. 

Study of a bundle of hexagonally ordered parallel 
F-actins with the charge of counterions distributed as 
columns of charge domains inside the bundle similar to 
the CDW picture of Ref. [1] showed that the configura- 
tion of minimum energy corresponds to helical sym- 
metry of F-actins with columns of charged regions in the 
middle of the separation between axes of neighboring F- 
actins. In this arrangement, sdls of F-actins and the 
charge domains of opposite sign are closer to each other 
in comparison with the other possible arrangement in 
which the columns of charged domains are located on 
the axes of the channels formed by three adjacent F- 
actins. This calculation shows that if we arrange the 
counterions according to the scenario of the charge den- 
sity wave presented in Ref. Q, the minimum energy of 
the system in our model corresponds to jj symmetry, not 
11 . Considering F-actin torsional modulus, overtwisting 
it from its native helical symmetry to for example H 
symmetry (3.8° over- twist per monomer) costs an en- 
ergy of ~ IGfcsT per monomer. The other symmetry, 
Iy , which corresponds to the global minimum of energy 
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in Fig. [TOl differs only 1.2° per monomer from F-actin 
native structure. Energy cost of such a deviation (under- 
twist) is in the range of thermal fluctuations strength and 
under-twisting of F-actin to this symmetry seems to be 
possible. In our simulations however, we found that the 
system prefers the F-actins to have local deviations from 
their native symmetries instead of regularly changing of 
the helical symmetry. 

Some important features of real F-actin macro- 
molecules have not been considered in our model. De- 
spite our assumption, it is known that G-actin monomers 
charge is distributed as discrete point charges on the sur- 
face of its domains. We assumed however in our model 
all the net charge of G-actin as a point charge in the 
center of its sdl. One should note that other models 
containing more details of F-actin may be obtained from 



simulations of F-actin in the atomistic scale [1^, . The 
other important feature of real F-actin that has not been 
considered in our model is the low value of its dielectric 
constant, ~ 2, relative to that of water, ~ 80. This differ- 
ence is known that can affect distribution of counterions 
in the vicinity of such macromolecules [47l | . 
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